1. Field of the Invention
The present invention relates to a method of computation and processing of the image in conical projection, for example in the manner of x-rays, of a three-dimensional (3D) object sampled in volume elements, as well as to a method of three-dimensional reconstruction of an object under study which makes use of this method of computation. The invention finds an application more particularly in an medical field in which the objects to be studied are patients' bodies subjected to a radiological examination. To this end, the invention is essentially concerned with three-dimensional reconstruction. However, it would also be applicable to methods of two-dimensional reconstruction. Similarly, the invention can be employed for visual display of a sampled volume. The object of the invention is to contribute to the production of representative images of 3D objects under study which are more sharply defined and also obtained more rapidly.
2. Description of the Prior Art
It is already known to make two-dimensional (2D) reconstructions of images of cross-sections of objects after (1D) acquisitions in radiological tomodensitometry performed in the cross-section to be imaged. The different generations of tomodensitometers have led to the use, in a third generation, of x-ray fan-beam tomodensitometers in which an x-ray point source illuminates in x-radiation a so-called multidetector comprising a plurality of cells aligned in the plane of the fan beam. The object to be studied is interposed between the x-ray source and the multidetector. Acquisition involves a series of illuminations. From one illumination to another, the assembly consisting of x-ray source and multidetector is rotated about the body to be studied. If s designates the longitudinal coordinate of one cell of the multidetector on said multidetector and if .theta. designates the angle of location of the x-ray source-multidetector assembly at the moment of illumination, there is accordingly obtained a series of radiological absorption measurements denoted by P(.theta.,s). If one designates as x and y the coordinates of a volume element of the object to be studied, in the cross-section concerned, and if one gives the notation f(x, y) to the linear attenuation function of the x-radiation which passes through the object, one may write: ##EQU1## In this expression, X designates the x-ray which arrives on the cell s when the x-ray tube - multidetector assembly is in the orientation .theta.. I.sub.0 is the incident intensity. This formulation is valid only in cases in which the source is monoenergetic; t is the curvilinear abscissa along the x-ray.
The application of the Radon theories had entailed the need to find the one-dimensional Fourier transform of P(.theta.,s) denoted by .rho..sub..theta. (u). By then finding the two-dimensional Fourier transform of f(x,y) denoted by .tau.(v,w), one could be led to identify .rho..sub..theta. (u) and .tau.(v,w) evaluated on the straight line of angle .theta.. It has been deduced from this that, by means of a two-dimensional Fourier transform reverse to .tau.(v,w) and obtained from all the transforms .rho..sub..theta. (u) (after changing of variables and interpolation), it was possible to compute the distribution of f(x,y) in the cross-section under study from the measurements P(.theta.,s).
In practice, one avoids Fourier transforms which lead to excessively long calculations and one makes use of a so-called filtering technique followed by back-projection. Filtering consists in computing the convolution product of the function P(.theta.,s) which is representative of the measurements, by means of a filter function q(s). This product is as follows: ##EQU2##
In this expression, t is a dummy variable of integration. The choice of the filter (which virtually consists in solving the same problems as those encountered at the time of changes of variables and of interpolation contemplated earlier) governs the quality of the images obtained. In the back-projection operation, the convolved value is assigned to all the volume elements (voxels) of the cross-section of the object under study which is located on the path of an x-ray (.theta.,s) concerned. In order to deduce the image therefrom, the different convolved values assigned are cumulated for each voxel.
Up to the present time, the acquisition of data and the 3D reconstruction of the structures examined has been carried out by displacing in translational motion the x-ray source - multidetector assembly along the body to be studied and by thus acquiring a plurality of adjacent 2D cross-section images in said body. However, this technique is wholly inadvisable in angiography in which a contrast product is injected in order to obtain contrast in the blood vessels. Injection of this contrast product has a traumatizing effect on the patient, particularly when it is repeated too often. It is for example prohibited to perform 256 injections of contrast product if it is desired to acquire images of 256 adjacent cross-sections in the body. Apart from this prohibition in angiography images, it must be acknowledged that the technique of 2D acquisition and corresponding 2D reconstruction is too time-consuming to be employed for reconstructing 3D objects. In fact, with a mean resolution of 256 points per 256 points in a cross-section image, the acquisition time with a tomodensitometer of current design is of the order of four seconds. The acquisitions which are necessary for 256 cross-sections would accordingly result in an examination time of nearly one half hour. This length of time is much too long to be endured by patients and for the public health system as a whole (general cost).
In theory, it is possible to generalize the Radon method by carrying out 3D acquisitions and by directly performing 3D reconstruction of the objects to be studied. By 3D reconstruction is meant the computation of a numerical volume in which memory cells placed at representative addresses of the voxels of the volume under study contain an item of information corresponding to the distribution of the (radiological) phenomenon being studied in the object. However, the Radon theory makes it necessary to acquire measurements corresponding to integration of the physical characteristics to be imaged in a sub-space assembly, or so-called hyperplane, the dimension of which is smaller by one unit than the dimension of the space under study which is to be reconstructed. In other words, in the case of a 3D space, it is necessary to have measurement results integrated on 2D hyperplanes. In point of fact, radiological acquisition with a point detector can only be an integration of the characteristics on a (1D) straight line: the x-ray. In the event that in practice, the trajectory of the source fails to satisfy the conditions laid down by theory, it will not be possible simply on the basis of a knowledge of the projections along all the straight lines (x-rays) which have served for acquisition, to compute the projections along all the hyperplanes which are necessary for reconstruction. In other words, points will be missing in the Radon space of the measurements. This latter will not be filled in a uniform manner. Artifacts will consequently appear in the resulting 3D reconstructions. Although two-dimensional multidetectors may at present be contemplated (arrangement of silicon photodetectors, use of radiological image intensifier screen), 3D reconstructions by this method must still be relegated to the rank of desirable objectives in view of the imperfection of results to which said method leads when one cannot supply a large number of data, as is the case in vascular applications.
Consideration has been given to a radically different algebraic reconstruction technique based on the following principles. The set of measurements P(.theta.,.tau.,s) acquired with a plane or concave 2-dimensional multidetector is already known. It is also known that there exists a continuous function f(x,y,z) which is representative of the radiological absorption phenomenon which it is endeavored to represent. It is sought to determine f by reconstruction. In practice, taking into account the fact that all the calculations are performed by treatments of the data-processing type, the knowledge of f at the exit is a sampled knowledge. The new approach has consisted in estimating f by means of a discrete function denoted by f and established a priori. For example, f can consist at the outset of a numerical volume in which all the voxels are established at one (or at zero). This discrete function is then projected on the multidetector as if the volume under study corresponded exactly to this discrete function. There is thus obtained an estimation designated as P.sub.i (f). This may also be written EQU P.sub.i (f)=H.sub.i f 3
In these expressions, i relates to a cell number i of the multidetector and H.sub.i represents the line i of the projection matrix H which corresponds to the cell i. The projection matrix H is independent of the body under study. Said matrix is dependent only on the projection geometry and constitutes the modelization of this latter. One then compares P.sub.i (f) (estimated) with P.sub.i (f) (measured). If there is no identity, and it is apparent that there is no identity at the outset, then f is corrected. In order to correct, there is performed on the object to be reconstructed a back-projection by the value of the difference between the measured projections and the computed projections. This is written EQU f.sub.k =f.sub.k-1 +.lambda..sub.k (H.sub.i *.epsilon..sub.i.sup.k)/.vertline.H.sub.i .vertline..sup.2 4
where H.sub.i is the adjoint operator of H.sub.i * and where .epsilon..sub.i.sup.k is the difference between the measured and computed projections at the iteration k-1.
This set of projections is repeated until the identity of P.sub.i (f) and of P.sub.i (f) is sufficient. This technique has been described in "Images reconstruction from projection" by G. T. Herman, Academic Press, 1980.
Moreover, if the projection matrix H is a binary matrix with 1s and 0s, the results are poor. Thus a first solution has consisted in approaching P by the curvilinear integral of f taken along a straight line D which connects the source S and the center of each cell on the multidetector. This amounts to considering, however, that each cell is of infinitely small size since, in this manner, neither the real surface of the cell nor the conical character of the projection on said cell is taken into account. Computation of this curvilinear integral may be carried out in several ways. The simplest method consists in approaching the curvilinear integral by a weighted sum of the values of the voxels traversed by the x-ray. There can thus be written: ##EQU3##
In this formula f.sub.j represents the absorption function in each of the voxels (the sampled volume elements) traversed, and h.sub.ij represents the weighting assigned to the value f.sub.j =f(x.sub.j,y.sub.j,z.sub.j) in the case of said voxel when it is traversed by a ray which arrives on the cell i. This weighting can be approached by the length of intersection between the voxel traversed and the x-ray. In this approximation, or so-called method of square pixels as described in the first document cited, the results produced are unfortunately of insufficient quality for the applications which are contemplated.
Another weighting computation has been conceived and is described in "An improved algorithm for reprojecting rays through pixel images" by P. M. Joseph, IEEE-MI, Volume 1, No. 3, Pages 192-196. The basic idea of this other computation is to approach the curvilinear integral of f by a sum relating to the interpolated values of f taken at the intersections of the straight line followed by the x-ray and of horizontal (or vertical) straight lines passing through the centers of the voxels visited or grazed by the x-ray which arrives at the center of the multidetector cell. The values of f at these points of intersection are estimated by means of a linear interpolation between two samples which are located in closest proximity, in the horizontal direction (or in the vertical direction) along the slope of the straight line of the x-ray.
A generalization of the method of square pixels has been proposed by K. M. Hanson and G. W. Wecksung, "Local basis - functions approach to computed tomography", Applied Optics, Vol. 24, No. 23, December 1985, Pages 4028-4039. In this generalization, a continuous function f(x,y,z) is constructed from samples f(x.sub.j,y.sub.j,z.sub.j) and the projection of this continuous function on the cell i of the multidetector is computed. In this generalization, the function f is defined as a linear combination of a set of basis functions b.sub.j (x,y,z) and we may write: ##EQU4##
In this formula, a.sub.j represents in the final analysis the description of the function f which is sought, the functions b.sub.j (x,y,z) being known. One usually makes use of basis functions which are local, repetitive and separable. "They are local" means that the support of each function b.sub.j is small compared with that of f. "They are repetitive" means that each function b.sub.j is deduced by translation from a single function b(x,y,z). This deduction is of the form EQU b.sub.j (x,y,z)=b(x-x.sub.j,y-y.sub.j,z-z.sub.j) 7
"They are separable" means finally that there exist three functions b.sub.x (x), b.sub.y (y), and b.sub.z (z) such that: EQU b (x,y,z)=b.sub.x (x).multidot.b.sub.y (y).multidot.b.sub.z (z)8
It is also worthy of note that, in the event that b.sub.x is equal to b.sub.y and to b.sub.z and that these functions are equal to a square-wave function, the so-called method of square pixels is again encountered. It may be stated that the generalization obtained by the method of basis functions makes it possible to revert to the other methods mentioned earlier but by giving them a more explicit mathematical foundation, with the result that the operations performed can be more readily understood. In practice, the best results are obtained with basis functions which are B-splines, cardinal sines, Hamming functions, or preferably Gaussian functions. With such basis functions, the results obtained are of better quality than those produced by the Joseph method, although of course at the cost of a longer computation time.
The operation which consists in projecting the estimation f of f actually consists in computing ##EQU5##
In this expression, g.sub.i designates a weighting function having a zero value (or a very small value) outside the angular sector formed by the source s and the cell i. The functions g.sub.i are representative of the conical support for illumination of the cell i. We may in fact write EQU h.sub.ij =.intg..intg..intg.g.sub.i (x,y,z).multidot.b.sub.j (x,y,z) dx dy dz 10
and the computation of p is then rcduced to calculation of all the weights h.sub.ij. Each weight h.sub.ij represents the contribution of the basis function b.sub.j to the projection of the object on the cell i.
Starting with a knowledge of samples f.sub.j of a function f of three variables, the aim of the invention is to permit numerical computation of the projection of f on a detecting cell i while ensuring both good quality of the result obtained and very short computation times. These two criteria (quality and speed) are essential in the reconstruction of a function f(x,y,z) from a set of 2D conical projections by means of algebraic techniques of reconstruction. Whereas it was not possible to ensure the desired quality and speed by means of the methods of projection of the prior art, the present invention makes it possible to do so.